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A dielectric solvation model is applied to the prediction of the equilibrium ioniza- 
, tion of liquid water over a wide range of density and temperature with the objective 

QQ I of calibrating that model for the study of ionization in water of organic acids, e.g., 

■ proteins and nucleic acids. The model includes an approximate description of the 

i polarizability of the dissociating water molecule. The calculated pK^^ are very 

sensitive to the value of the radii that parameterize the model. The radii required 
for the spherical molecular volumes of the water molecule in order to fit the exper- 
imental ion product are presented and discussed. These radii are larger than those 
commonly used. They decrease with increasing density as would be guessed but 
1^ . the rate of decrease is slight. They increase with increasing temperature, a varia- 

I tion opposite to what would be guessed if radii were strictly viewed as a distance 

of closest approach. The molecular theoretical principles that might provide an 
explanation of the thermodynamic state dependence of these radii are discussed. 



1. Introduction 

The chemical reaction H2O = H0~+ H"*" is of central importance to a 
variety of practical problems including those associated with chemical pro- 
cesses in aggressive aqueous environments such as supercritical water [1, 2]. 
It has been considered the most important chemical reaction in aqueous so- 
lution [3]. 

The equilibrium point of this process is often discussed on the basis 
of the ion product Kw = [HO~][H"'"]. Measurements of Kw for water at 
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densities below about twice the critical density are not available though 
convincing model estimates of this property have been given for the higher 
temperatures and pressures characteristic of supercritical water processes 
[4, 5]. Moderately high temperatures are required in order to examine Kw 
over a continuous and extended domain of thermodynamic states, including 
densities lower than those liquid densities most commonly encountered. 

This report gives theoretical calculations of Kw over a wide range of 
thermodynamic states. The calculations here employ a dielectric solvation 
model similar to that used by Pitzer and collaborators [4, 5]. As an 'equation 
of state' model the present results are not expected to improve upon that 
prior work. But the present contribution is helpful for several reasons that 
focus on the basis of the theory. The model used here is simple and the 
connection with molecular theory is simple too [6-10]. Thus the results here 
will facilitate direct molecular-scale verification of a number of intermediate 
quantities. This topic is of recent theoretical interest [11-14]; we expect 
comparable molecular-scale theory and results to be available in the near 
future. 

Since the model we exploit is common in studies of the charge states of 
protein and nucleic acid molecules in aqueous solutions, the molecular scale 
verifications are expected to improve our understanding of those biophysical 
problems too [15-23]. 

A concern with dielectric solvation models is the assignment of cavity 
radii or volumes. Cavity volumes should depend on the thermodynamic 
state [8, 10, 24, 25]. This is important to the extent that derivatives of the 
solvation free energy with respect to thermodynamic state, i.e. temperature, 
pressure, and composition, are required. The molecular theoretical principles 
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upon which those radii might be determined have been identified [7, 8, 10]. 
However, this issue can also be studied empirically by analyzing the solvation 
free energy computed on the basis of the dielectric solvation model over 
a sufficiently wide range of thermodynamic conditions. The data on the 
chemical equilibrium that is studied in this report are more extensive than 
for most reactions in aqueous solutions. The water self- ionization reaction 
is thus suitable for the principal motivation of this work: the delineation of 
the region where the model works satisfactorily and a calibration of the radii 
required to describe ionization of organic acids. 

The next section gives technical details of the theoretical issues and 
methods. We then present the results and a discussion. 

2. Methods 

The dielectric solvation model used here is formulated from the following 
physical ingredients. Attention is focused on a solute of interest and a solute 
volume is defined. For liquid water under the most common conditions it 
is known that the van der Waals volume of the molecule is a reasonable 
initial choice for the molecular volume. Partial charges describing the solute 
electric charge distribution are positioned within this volume. The solvent 
is idealized as a continuous dielectric material with the measured dielectric 
constant s. The solvent is considered to be excluded from the solute volume 
and that molecular volume is assigned an internal dielectric constant Sm', 
when £^ > 1 the model thereby includes an approximate representation of 
the molecular polarizability [26] . In all the calculations below we have used 
Sm ~ 2.42 which corresponds approximately to the known polarizability of 
the water molecule (1.44 A^) when a spherical molecular volume of radius 
1.65 A is assumed. Note, however, that when adjustments of the radii are 
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considered below the value for Sm is not adjusted. 
The equation to be solved for the model is 

V«£(r)V$(r) = -47rp/(r), (1) 

where Pf{r) is the density of electric charge of the solute molecule, the func- 
tion £(r) gives the local value of the dielectric constant, and the solution 
$(r) is the electric potential. To solve this equation, we first cast it as an 
integral equation, e.g., 

$(r) = $(°)(r)/e^ + J G(°)(r,r') • V'$(r')rfV. (2) 

V 

Here G^'^\r, r') is the Green function for the Poisson equation with £(r) = 1 
and $('^^(r) is the electrostatic potential for that case. This equation is 
correct both for a localized distribution p/(r) and zero boundary data on 
a surface everywhere distant and for periodic boundary conditions on a cell 
of volume V. G^^\r,r') is different in those two cases as is $W(r). For 
the unbounded cases treated here G^^\r,r') = l/|r — r'|. This equation is 
not the only such form that can be solved but it is the basis for the most 
popular boundary element approaches for numerical solution of this physical 
model [27-33] . Because of this we defer a more detail discussion of integral 
formulations of Eq. 1. 

The integrand of Eq. 2 is concentrated on the interface between the so- 
lute volume and the solvent. The interpretation of the integral of Eq. 2 is 
that it superposes contributions to the electrostatic potential due to charge 
induced on the surface. We can then use boundary element ideas to solve it. 
We have used sampling methods, here principally those based upon quasi- 
random number series [34-36], to evaluate the surface integral rather than 
more specialized methods. Advantages of such approaches are that they 
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are simple, yet permit systematic studies of numerical convergence and ex- 
ploitation of systematic coarse-graining. The geometries considered here are 
particularly simple, so these applications are not a stringent test of numerical 
methods. But an elegant statement of the advantages of boundary element 
approaches in the this context can be found in a paper of Yoon and Lenhoff 



In this application we solved Eq. 2 for the molecule in the liquid with the 
experimental dielectric constant and in isolation with the external dielectric 
constant assigned a value of e = 1. The solutions for these cases are denoted 
by ^i{r) and $i,(r), respectively. Remember that the 'internal' dielectric 
constant parameter is fixed at Em — 2.42 throughout. Then we construct 



the standard chemical potential difi'erence used below. This difference in- 
volves only the indirect contributions of Eq. 2, i.e. an infinite self-energy 
is eliminated by this subtraction. Here p/(r) is a sum of partial charges. 
Therefore the integral in Eq. 3 is a sum over those partial charges. 

With applications of such an ad hoc model, some arbitrariness is practi- 
cally unavoidable. For these calculations the most important features of the 
method that are arbitrary are two: (i) the chemical equilibrium and chemical 
species that are treated; and (ii) the intramolecular volumes assigned. 

The chemical reaction that we study in order to calculate Kw is 



[30]. 



A/i = 




(3) 



V 



2 H2O = HsO+ + HO-. 



(4) 



This arbitrary choice is made on the basis of the physical guess that obser- 
vation of a 'free' proton would be less likely than the case where an excess 
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proton is associated with an oxygen atom in roughly the same circumstances 
that an H-0 atom pair are configured in the water molecule. Atoms more dis- 
tant from an 'excess' proton than the nearest oxygen should be of secondary 
importance. We then consider the equilibrium ratio 

^^[HsO-nHO-]_ (5) 

We express this quantity as 

K = (T) X exp {- {Af^Ho- + ^f^HsO+ - '^AfiH.o) /RT} , (6) 

where K^^^ (T) is the equilibrium constant in the absense of intermolecular 
interactions, i.e. for the ideal gas [11, 12]. Accounting for the mass balance 
of the Eq. 4 we then find 



Kw = [H30+][H0-]=p'K 



l-AK 



(7) 



where p is the formal density of the water. 

The ideal gas contribution K^^^ (T) can be obtained from standard for- 
mulae and information on the molecular structures of the species considered. 
The information used here is collected in Table I [37-39] . 

What remains is the evaluation of the excess chemical potentials ap- 
pearing in Eq. 6 for the equilibrium ratio K. This brings us to the second 
arbitrary feature of our calculation, the assignment of molecular volumes. 
Here we assumed that each molecular volume is a sphere centered on the 
oxygen atom and that the radii of the spheres for all species were equal. 
This constraint is not necessary. However, some fitting to experimental re- 
sults will be required and limiting the number of fitting parameters seems 
prudent. It might be hoped that the fitted radii would validly describe ioniza- 
tion of organic acids too; if the radii of protonated and deprotonated oxygen 
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atoms were not transferrable to different but similar chemical environments 
with acceptable accuracy much of the usefulness of the model would be lost. 
In addition, the molecular underpinnings of the dielectric solvation model 
have been clarified over recent years; for solute species with spherical short- 
ranged (non-electrostatic) interactions with the solvent, we can hope that 
sufficient molecular information will become available to make applications 
such as this more geniunely predictive [7, 8, 10]. 

We note the precedent for cations to be assigned slightly larger radii 
than the isoelectronic anions [40, 41]; several explanations of this have been 
offered [8, 40, 41]. Here we assume that such complications are sufficiently 
slight that we can ignore them. 

Finally we must specify the charge distribution p/(r) of the chemical 
species treated and the dielectric constant. For the former quantity we used 
atom centered partial charges following the method of Breneman and Wiberg 
[42] . The dielectric constant utlized was that given by the empirical equation 
of Quist and MarshaU [43]. 

3. Results and Discussion 

Our results for Kw are shown in Fig. 1. The densities treated range over 
a factor of three including the triple point density. The lowest density shown 
is approximately 0.4 g/cm^. Temperatures encompass all of the traditional 
liquid range from the triple temperature to above the critical temperature. 
Fig. 2 shows the empirical radii that resulted from fitting the model result 
to experimental values given by Marshall and Franck [3]. The magnitudes 
of these radii are reasonable but they are significantly larger than would 
be most commonly guessed [41]. For the upper part of this density range, 
the radii are decreasing functions of density at constant temperature. But 
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the variations of the radii with density at fixed temperatures are shght; the 
decreases in the radii are 3-7 times smaller than the rule SR/R ~ —dp/3p 
based upon the Onsager estimate of the radius, AnR^p/S = 1 [24, 44]. The 
variation observed for R with temperature at fixed density is not accounted 
for by the Onsager radius. Moreover, a natural guess for the variation of R 
with temperature at constant density is not so clear; if R were identified as 
a distance of closest approach of the solvent molecules to the solute then R 
would be expected to decrease with increasing temperature at fixed density. 

The sensititivity to the parameter R is an important issue for applica- 
tions of the model. In the present case, pK^i^ can change by 1 unit if the value 
of R is changed by 1% from the value that fits the experiment. A change of 
2% in R from the optimum value can lead to a change of more than 3 units in 
pKw Computational searches for a set of three radii [Rn-iOi RhzO+i Rho-) 
that could satisfactorily fit the whole behavior in Fig. 1 without thermo- 
dynamic state dependence were unsuccessful. For example consider the 
isotherm at T = C of Fig. 1 for which the total change in pKvi^ is about 
1.7. The set of three radii that we found to give the best fit to those data still 
lead to errors of -1-0.57 and -0.55 at the beginning and the end, respectively, 
of the isotherm shown. The set of three radii that produced the best fit to 
the data over the whole of Fig. 1 lead to typical errors in pK^^ greater than 
4 and those radii were not chemically reasonable. 

How the model corresponds most naturally to molecular theory is known 
[7, 8, 10] and this correspondence provides an avenue for theoretical de- 
termination of R. The corresponding molecular theory is the second-order 
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cumulant approximation 



A// ^ Alio + ( "^(j) 




Here <^(j) is the electrostatic interaction potential energy coupling the solute 
to solvent molecule j. The brackets (• • •)q indicate the thermal average in 
the absence of those electrostatic couplings and A/iq is the excess chemical 
potential of the solute at infinite dilution again in the absence of electro- 
static interactions. When the solute is an infinitely dilute second component 
this molecular approximation is perturbation theory through second order 
in the electrostatic interactions. When the molecule is not literally an in- 
finitely dilute second component the 'infinite dilution' restriction means that 
one molecule is distinguished for the purposes of calculation. This is still 
a natural theory but the medium now contains a non-zero concentration 
of molecules mechanically identical to the 'solute.' The medium properties 
non-perturbatively reflect that fact. From the perspective of this molecular 
theory, the dielectric solvation model application neglects the zeroth and first 
order terms, makes an estimate of the second-order term, and neglects all 
succeeding contributions. 

Careful calculation on a molecular basis of the terms in Eq. 8 would be 
a helpful next step for clarification of these theoretical models. Jayaram, et 
al. [45] and Rick and Berne [46] have tested dielectric solvation models by 
computer experiment. Their results provide, in principle, the information 
requested here but they were not analyzed from the point of view of the 
present goals for this system. Those results together suggest that hydrogen- 
bonding interactions are border-line cases for dielectric solvation models. 
Although the present application and its motivation in ionization of organic 
acids is of wide interest, it is likely a severe test for these models. Thus 
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the empirical radius R determined here probably subsumes other errors that 
would have been encountered had the molecular theory been implemented. 
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Figure Captions 

Figure 1: The water ion product ([mol/1]^) as a function of density p (g/cm^) 
and temperature over a range of liquid thermodynamic states. The re- 
gion underneath the domed-shaped curve is the hquid- vapor coexistence 
region. 

Figure 2: Empirical radii R, as a function of density and temperature, required 
for the dielectric solvation model to fit the experimental results. 
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Table I. Properties used in the calculation of molecular partition 
functions 
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5796.13 


2736.17 
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4267.23 
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4304.29 




Qv6 




4315.57 






13.6745 


8.3841 


33.9547 


Qb 


20.9686 


13.4501 




Qc 


39.1187 


13.9630 




E 

^9 


-1032.98 


-1203.55 


-643.127 



Geometry from Ref. 38. 
^ Geometry from Ref. 39. 
Rotational symmetry number. 

Vibrational temperatures as defined in Ref. 37, in kelvin. These prop- 
erties were calculated using SCF approximation, a 6-31g* basis set, and 
the geometries of Refs. 38 and 39. 

" Rotational temperatures as defined in Ref. 37, in kelvin. These prop- 
erties were calculated in the same way as the vibrational temperatures 
were. 

Electronic ground state energies from Ref. 39 relative to the energies of 
the separated isolate atoms, in kcal/mol. 



15 




-0.4 -0.3 -0.2 -0.1 0.0 0.1 

Logio(p) 




p (g/cm3) 



